Integrating growth and survival models for flexible estimation of size‐dependent survival in a cryptic, endangered snake

Abstract Estimates of demographic rates for animal populations and individuals have many applications for ecological and conservation research. In many animals, survival is size‐dependent, but estimating the form of the size–survival relationship presents challenges. For elusive species with low recapture rates, individuals’ size will be unknown at many points in time. Integrating growth and capture–mark–recapture models in a Bayesian framework empowers researchers to impute missing size data, with uncertainty, and include size as a covariate of survival, capture probability, and presence on‐site. If there is no theoretical expectation for the shape of the size–survival relationship, spline functions can allow for fitting flexible, data‐driven estimates. We use long‐term capture–mark–recapture data from the endangered San Francisco gartersnake (Thamnophis sirtalis tetrataenia) to fit an integrated growth–survival model. Growth models showed that females reach longer asymptotic lengths than males and that the magnitude of sexual size dimorphism differed among populations. The capture probability and availability of San Francisco gartersnakes for capture increased with snout–vent length. The survival rate of female snakes exhibits a nonlinear relationship with snout–vent length (SVL), with survival flat between 300 mm and 550 mm SVL before decreasing for females between 550 mm and 700 mm SVL. For male snakes, survival decreased for adult males >550 mm SVL. The survival rates of the smallest and largest San Francisco gartersnakes were highly uncertain because recapture rates were very low for these sizes. By integrating growth and survival models and using penalized splines, we found support for size‐dependent survival in San Francisco gartersnakes. Our results have applications for devising management activities for this endangered subspecies, and our methods could be applied broadly to the study of size‐dependent demography among animals.

Vital rates can vary among individuals within a population depending on characteristics such as age (Gamelon et al., 2016;Zheng et al., 2007), sex (Ergon & Gardner, 2014), and size (Hyslop et al., 2012). Identifying how vital rates are influenced by individual characteristics has value for answering both basic and applied ecological questions. Interindividual variation in vital rates can reveal how selection acts to favor certain traits (Letcher & Horton, 2008) or differs among sexes (Schulte-Hostedde et al., 2002;Toïgo & Gaillard, 2003). In wildlife conservation, quantifying variation in vital rates among individuals is valuable for designing management actions that will have the greatest effect on population growth (de Kroon et al., 2000). One trait that can explain variation in vital rates for many animals is individual body size (Sauer & Slade, 1987a).
Body size is a fundamental trait that influences metabolism (Kleiber, 1947), lifespan (Speakman, 2005), and age-at-maturity (Stearns, 1984). Growth (Armstrong & Brooks, 2013;King et al., 2016), survival Sauer & Slade, 1987b), and fecundity (Rose, Ersan et al., 2018;Weatherhead et al., 1995) can all be size-dependent in animals. Researchers have developed a variety of models of individual growth (Andrews, 1982;Armstrong & Brooks, 2013;Keevil et al., 2021) and survival (Lebreton et al., 1992) to estimate these vital rates and explore how covariates explain individual, spatial, and temporal variation. Integrating growth and survival models can reveal insights into life-history patterns that would not be identifiable from either model alone (Reinke et al., 2020). For example, estimating the effect of individual size on survival probability in a capture-mark-recapture (CMR) model requires a value for the size of each individual during each sampling period. For elusive animals with low recapture probabilities, it is impossible to measure each individual during each sampling period. If individuals cannot be measured during each sampling period, a growth model can be used to estimate size when it is unobserved (Bonner & Schwarz, 2006;Schofield & Barker, 2011). Body size can also influence capture probability Schofield & Barker, 2011) and availability for capture (i.e., emigration; Riecke et al., 2018); therefore, estimates of the size-survival relationship could be confounded unless data are analyzed in a CMR model that separates the observation, emigration, and survival processes (Kendall et al., 1997;Koons et al., 2009). Integrating a model of individual growth into a CMR model using Bayesian methods enables recapture and survival probabilities to be modeled as a function of individual size, using estimated size (with uncertainty) when no empirical data are available (Reinke et al., 2020;.
Determining the shape of the functional relationship between vital rates and individual covariates presents another challenge.
Studies often assume a linear relationship between size and survival on the logit scale (Hansen et al., 2015;Wallace et al., 2013). There might be no reason to assume a priori a linear relationship between size and survival, and enforcing a linear function could lead to a misspecified model and researchers concluding any relationship is weak or nonexistent. For example, survival might increase with size up to a point before approaching an asymptote or declining for individuals above a certain size threshold (Brown & Weatherhead, 1999;Doak et al., 1994). If a linear relationship is not the only plausible option, one approach is to fit and compare multiple parametric functions (e.g., linear, quadratic, asymptotic) in a model selection framework . This does not completely resolve the issue, because each parametric function enforces a particular shape on the size-survival curve, which can lead to unrealistic behavior in the tail, such as a symmetrical decrease in survival for large and small individuals with a quadratic function (Gimenez et al., 2006). Choosing a suboptimal parametric function when estimating relationships between individual covariates and vital rates can have cascading effects when vital rate functions are used to construct a demographic population model (Dahlgren et al., 2011), which could hinder such a model's ability to provide insight into management of sensitive species. Nonparametric spline functions present an appealing option when the shape of the functional relationship is unknown, because the data dictate the shape of a flexible curve (Gimenez et al., 2006).
To prevent overfitting, penalized splines draw coefficients from a shared distribution, shrinking spline parameters toward zero if there is little support for their effect (Crainiceanu et al., 2005). Fitting splines can be thought of as a form of data exploration (Tredennick et al., 2021), rather than strict hypothesis testing about how vital rates vary with individual traits. In practice, spline regressions have been used to derive smooth functional relationships between vital rates and covariates in demographic models (Robinson & Wilson, 2021;Schwinn et al., 2017).
Using nonparametric splines in conjunction with CMR models has the potential to reveal valuable information on the demography of species of conservation concern. Here, we fit an integrated growth-survival model to long-term CMR data from an endangered snake subspecies, the San Francisco gartersnake (Thamnophis sirtalis tetrataenia). Capture probabilities for this species are low (Halstead et al., 2011;Reeder et al., 2015); therefore, fitting a growth model was necessary to assess how individual size and sex affected survival, recapture probability, and availability for capture. We also evaluated whether and how individual  2 | ME THODS

| Study species and sites
The San Francisco gartersnake ( Figure 1) is a subspecies of the common gartersnake (T. sirtalis) that occurs in San Mateo and northwestern Santa Cruz counties in California, USA, and is listed as endangered under the federal Endangered Species Act (ESA, 1973;U.S. Fish & Wildlife Service, 1967) and by the state of California (Fish & Game Commission, 1971). Common gartersnakes occur in various types of habitat (Fitch, 1941), but the San Francisco gartersnake subspecies is reliant on wetlands, riparian forest, and adjacent upland habitat (Barry, 1994;U.S. Fish & Wildlife Service, 2006). San Francisco gartersnakes primarily feed on amphibians (Kim et al., 2021) including native Sierran treefrogs (Pseudacris sierra), California red-legged frogs (Rana draytonii), and Pacific newts (Taricha spp.).
San Francisco gartersnake populations have become more isolated as urbanization in the species' range has altered wetlands and surrounding upland habitats (U.S. Fish & Wildlife Service, 2020;Wood et al., 2020). The minimum reported size at sexual maturity for female San Francisco gartersnakes at one population (not studied here) was 368 mm snout-vent length (SVL) with a modal SVL of 420 mm for gravid females (Reeder et al., 2015). Females are expected to mature at 2-3 years of age (U.S. Fish & Wildlife Service, 2020).
We sampled for San Francisco gartersnakes at a total of five sites ( Figure 2). We obscure the precise locations of these sites because illegal collection is a concern for San Francisco gartersnakes (U.S. Fish & Wildlife Service, 2006). We sampled two sites in San Mateo County, California, to collect CMR data for estimating survival (Site C and Site I). We sampled at Site C between 2007 and 2020 for an average of 51 ± 11.3 (SD) (range = 21-64) days (i.e., secondary occasions) each year and at Site I between 2014 and 2017 for 53.5 ± 7.5 (48-64) days each year. In addition to our two long-term study sites, we sampled three additional sites (Sites N, P, and S) for San Francisco gartersnakes in 2018 and 2019 for 37.5 ± 5.5 (30-45) days each year, as part of a study on population abundance and genetic diversity (Wood et al., 2020). Sites N, P, and S are all protected from human development; sites N and P are located in coastal parks and have public access, whereas Site S is located 10 km inland and is closed to the public. We use the data from these three additional sites for modeling individual growth. Despite Site S comprising two separate clusters of trap arrays, we treat Site S as a single site based on the close proximity between our trap arrays (7.7 km), the lack of a geographical barrier to movement, and snakes from this area belonging to the same regional genetic cluster (Wood et al., 2020).

Site C was an open-space preserve in coastal southern San Mateo
County (mean elevation: 110 m, range = 50-190 m) that was off-limits for recreation but subject to cattle grazing during the 20th century until 1998, and again from 2015 to 2020. Site C is regularly covered in coastal fog (Torregrosa et al., 2016), and the diurnal temperature ranged between 3.6 and 33°C (mean = 14.9°C, SD = 3.5) during our sampling.

| Data collection
Most sampling for San Francisco gartersnakes took place in April and May, although occasional surveys took place in March and later in the summer in some years (Table S1). We sampled all sites for San Francisco gartersnakes using trap arrays with 8-to 12-m-long drift fences made from tempered hardwood boards that were buried 5-7 cm in the soil and stood 40 cm high. Each drift fence had a total of four wooden box funnel traps (measuring 30 cm × 40 cm × 23 cm), with two traps placed flush against either side of the fence at each end. We installed trap arrays in upland habitat adjacent to wetlands; all trap arrays were <215 m from wetlands (Table S2). Sites I, N, S, and P were each sampled with 48 funnel traps deployed along 12 drift fences, and Site C was trapped with 96 funnel traps along 24 drift fences to cover the larger area of the site (Table 1). Trap arrays were installed at the same location each year, with minor changes if the wetland area expanded to inundate the normal location of a trap array. In addition to funnel traps, artificial cover objects were deployed at sites C and I, and snakes were captured by hand opportunistically during trap checking at all sites. The number of cover objects active each year at Site C was dependent on environmental conditions and space use of the cattle (Table S1). Funnel traps were checked daily while active, and we measured snout-vent length (SVL) to the nearest millimeter (by laying snakes on a meter stick) and mass to the nearest gram (using Pesola ® spring scales) of all captured snakes. We determined snake sex by cloacal probing and marked snakes by branding ventral scales with a unique code (Winne et al., 2006). We measured air temperature in the field while checking traps, and recorded whether precipitation fell in the prior 24 h.
We acquired data on mean daily relative humidity from interpolated gridMET data (Abatzoglou, 2013). For some dates when air temperature was not collected in the field, we used gridMET data to impute the air temperature, based on a linear regression model relating air temperature recorded in the field and gridMET mean daily temperatures for that site (Appendix S1).

| Growth model
We fit a von Bertalanffy (VB) growth model based on Armstrong and Brooks (2013) to snake SVL data from recaptures. Unlike Armstrong and Brooks (2013), we did not incorporate biphasic growth, because preliminary model fitting indicated there was no evidence for a fixed change in growth rates as snakes increased beyond a particular size threshold. This growth model takes the general form (Equation 1) where EL i , t is the expected SVL of snake i at time t, a i is the asymptotic SVL of snake i, k i,t is the somatic growth rate of snake i at time t, and Δt is the number of days between captures, including the inactive season.
Although the exponential component of the VB growth model is usually parameterized in terms of k (Wang, 1998), we followed Armstrong and Brooks (2013) and substituted k/a for k, which allows for individual variation in growth rate to be independent from individual variation in asymptotic size. The observed SVL (L i,t ) is drawn from a normal distribution centered on EL i,t with random variation around the expected size at time t (ε t ) that encompasses both measurement error and individual variation around expected size (Equation 2).
We modeled log(k) as a function of sex (fixed effect) and random effects of individual (i), site (s), and year (t) (Equation 3).
The variable µ f is the mean value of k for a female snake, β m is the effect of maleness on k, m i is a binary indicator of whether an individual is male (1) or not (0), and α sex,s , ζ t , and ι i are site (by sex), year, and individual random effects, respectively.
We modeled asymptotic length (a i ) as a function of sex and site (Equation 4). The mean estimate of a for a female snake is υ f , θ m is the effect maleness on a, m i is the binary male indicator defined above, and η sex , s is the sex-and site-specific random effect on a.
We defined the site, year, and individual random effects with normal distributions with a mean of zero and standard deviation shared among members of that cluster level ( Table 2). The site random effects on k and a were each sex-specific, which allowed these growth parameters to vary among sites in different ways for females and males, because interpopulation differences in sexual size dimorphism have been found in other subspecies of T. sirtalis (Krause et al., 2003). We then used the estimated a and k parameters and Equation 1 to project size at age for San Francisco gartersnakes using a neonate size at birth of 165 mm SVL based on measurements reported by Barry (1994).
We implemented the VB growth model using Markov chain Monte Carlo (MCMC) sampling in JAGS (Plummer, 2003) accessed through R version 4.0.5 (R Core Team, 2021) using the "runjags" package (Denwood, 2016). We ran the VB growth model on four chains for 2,500,000 sampling iterations each, after discarding the initial 20,000 iterations as burn-in. The resulting output was thinned by a factor of 100, resulting in a final sample of 100,000 model iterations used for inference. We evaluated model convergence by inspecting MCMC trace plots and calculating the Brooks-Gelman-Rubin statistic (R; Brooks & Gelman, 1998). All parameters had R < 1.01, and chains showed no evidence of lack of convergence. We tested the goodness of fit of the VB growth model using a posterior predictive check, by comparing the sum of squares of residuals for observed snake SVLs at recaptures compared with expected values, and repeating this calculation of residuals for replicate snake SVLs generated from the model (Keevil et al., 2021). There was no evidence of lack of fit in the scatterplot of residual values ( Figure S1, Appendix S1), and the Bayesian p-value was .48, indicating adequate fit of the model to growth data (Kéry & Schaub, 2012). The connection between parameters in equations 1-4 above and JAGS growth model code is presented in Table S3.

| Robust-design model
We fit a hierarchical, multistate, robust-design model (Riecke et al., 2018) to our CMR data of San Francisco gartersnakes from two longterm study sites (Site I and Site C) to estimate the annual apparent survival (ϕ), probability of being available for capture (γ), and daily recapture probability (p) at each site, while accounting for variation due to individual size. Our robust-design model conditioned on initial capture and is a multistate version of a Cormack-Jolly-Seber model (Cormack, 1964;Jolly, 1965;Seber, 1965 were not measured during their first capture, their initial SVL was imputed based on a normal distribution with the mean and standard deviation of SVL for that sex. We set informative priors for the VB growth parameters in the robust-design CMR based on the posterior distributions from the VB growth model (Table 2), with slightly wider uncertainty.
Under the robust design, sampling consists of widely separated primary periods (t) each composed of multiple secondary periods (j) separated by small temporal intervals; populations are assumed closed between secondary periods (i.e., no mortality or migration), while populations are assumed to be open (i.e., mortality and migration can occur) during the larger intervals between primary periods (Pollock, 1982). We assumed that individual snakes are in one of three states during each primary period: available for capture, temporarily unavailable, or permanently unavailable (dead or permanent emigrant). We used a random emigration formulation of the robustdesign model that assumed availability for capture in one primary period did not depend on an individual's state in the previous primary period. We assumed random emigration because our sampling was focused near aquatic habitats where San Francisco gartersnakes are likely to forage for amphibian prey each spring (Kim et al., 2021), and therefore, we did not expect an individual's availability for capture in one year to influence their availability the next year. We created our robust-design survival model in a Bayesian framework following the Riecke et al. (2018), which allows for unbiased estimates of apparent survival (ϕ), availability (γ), and capture probability (p). The key innovation in the model of Riecke et al. (2018) is that data likelihoods are modified such that only individuals known to be available for capture in a given secondary period affect the capture probability for that secondary period.
We included potential effects of snake SVL on ϕ and p by fitting penalized splines (Eqs. 5,6). Given the difficulty of obtaining precise estimates of availability on-site for capture without auxiliary data (Bird et al., 2014), and based on preliminary model fitting in which a spline function for γ approximated a linear relationship with SVL, we fit a simple linear effect of SVL on availability (Equation 7). In equation 5, π ϕ,sex,s is the sex-and site-specific intercept of apparent survival on the logit scale, and year-specific random effect on ϕ is ξ t .
In equation 6, ο p,sex,s is the sex-and site-specific intercept of p on the logit scale, and ω s,t is the site and year random effect on p. In addition to the effect of snake size on p, we included potential linear effects of air temperature, recent precipitation (rain the last 24 h before checking traps), relative humidity, and day of year on p, where the covariate matrix is W, and the vector of coefficients is δ p (Equation 6). In equation 7, ς γsex,s is the sex-and site-specific intercept of availability on the logit scale, λ γ is the slope of the relationship between snake SVL and availability for capture, and s,t is the site and year random effect on γ. We constrained the effects of survey covariates (e.g., air temperature, day of year) on capture probability to be linear functions because our two-month sampling period did not cover the full seasonal variation in weather at our study sites. Although we believe that using linear functions was a reasonable assumption for this short time period, it is likely that some of these covariates might have more complex relationships with capture probability when evaluated over a wider range of values (e.g., capture probability might be highest for an optimal temperature range and decline temperatures above and below that range).
We fit low-rank thin plate splines following the methods of Crainiceanu et al. (2005). We divided the range of SVL values into G = 5 knots (where g indicates an individual knot) and fit spline functions (degree = 3) to each knot (Equation 8), where q is a fixed effect of SVL, b g are random coefficients, κ are fixed knots along the range of observed SVL values, κ g is the quantile of SVLs that matches g/G + 1, and (xκ g ) + represents max(0, xκ g ), and d is the degree of the spline. The splines are penalized because each b is drawn from a normal distribution (N) with a shared standard deviation parameter (σ b ), making each b parameter analogous to a random effect (Equation 9). The penalty that comes from drawing splines from a shared distribution shrinks the b parameters toward zero and limits the smooth function from assuming too complex a shape.
TA B L E 2 Parameters, prior distributions, and summary of posterior distributions from the von Bertalanffy growth model for San Francisco gartersnakes (Thamnophis sirtalis tetrataenia) captured at five sites from 2007 to 2020. The prior "N(mean, SD)" represents a normal distribution with mean and standard deviation (SD), "Exp(scale)" is an exponential distribution with a scale parameter, and "-" indicates a derived parameter. Mean and SD are summary statistics from the posterior distribution for each parameter, and 2.50% and 97.50% indicate percentiles of the posterior distribution For a given response variable (ϕ or p), we used a shared spline function, f(x) among sites and both sexes. In other words, although the intercept could differ among sexes and sites, we assumed that the size-survival function would be the same for each site-sex group, likewise for the size-capture probability function. We also assumed a shared slope among sites and sexes for the size-availability functions. We made the assumption of shared shape of size effects among sites because most of our CMR data came from Site C, and we lacked sufficient data to estimate independent functions for snakes from Site I. We tested for potential correlation among annual growth and survival rates by using a covariance matrix to define annual random effects for these parameters (Appendix S1). Because a small number of snakes were of unknown sex, we imputed the sex of these snakes using a Bernoulli prior for a binary indicator of male sex with a probability of 0.5. We assumed a fixed 1:1 sex ratio for snakes of unknown sex based on previous estimates of the sex ratio in our study populations (Wood et al., 2020).
We ran the robust-design model using Markov chain Monte Carlo (MCMC) sampling in JAGS (Plummer, 2003) accessed through R version 4.0.5 (R Core Team, 2021) using the "runjags" package (Denwood, 2016). We ran the robust-design model on 10 chains for 50,000 iterations per chain after discarding a burn-in of 10,000 iterations. We thinned the resulting chains by a factor of 10, resulting in a final MCMC sample of 50,000 iterations. As for the VB growth model, we evaluated model convergence by inspecting trace plots and calculating R ; all R were ≤1.05, and chains showed good mixing.
We also evaluated the goodness of fit of the robust-design model by simulating replicate capture data from the model and comparing how replicate data and real data deviated from expected values of the number of recaptures . We used the Freeman-Tukey statistic as a measure of the discrepancy between expected values and real and simulated data (Brooks et al., 2000).
We made a scatter plot of the deviation for replicate and real data and calculated a Bayesian p-value to determine whether the robustdesign model provided adequate fit to the data ( Figure S2). The Bayesian p-value for the robust-design model was 0.48, indicating adequate fit to the data. The connection between parameters in equations 5-9 above and JAGS CMR model code is presented in  Figure 3). There was weak support for male San Francisco gartersnakes approaching their asymptotic SVL faster than females on average (posterior probability [k m > k f ] = 0.80; Table 2). Estimates of k were similar between males and females for each site except Site C, where k m,C had 0.96 posterior probability of being greater than k f,C ( Figure S3). Individuals of both sexes grew at similar rates until reaching 2 years of age, at which point the growth trajectories diverged, with male SVL rapidly approaching its asymptote by ages 4-5 years, and female SVL continuing to increase before beginning to plateau around 7-8 years (Figure 3). There was more variation in male asymptotic SVL among sites than female asymptotic length; in particular, male asymptotic SVL was greater at Site I than at Site C ( Figure S4). Greater variance in k was explained by differences among sites and individuals than differences among years (Table 2).

F I G U R E 3
Growth trajectory for female (black line, gray shading) and male (red line, red shading) San Francisco gartersnakes (Thamophis sirtalis tetrataenia) as a function of age. Each trajectory is based on the mean parameter estimates from a von Bertalanffy growth model (i.e., the values for an "average" site). The solid line represents the median projected size at each age, and the shaded region is the 95% credible interval. Growth projections assume that neonates of both sexes have mean SVL = 165 mm (SD = 10 mm) based on measurements reported in Barry (1994)  year, 35 were captured in two years, and 6 were captured in 3 years.
Daily recapture probability (p) was positively related to air temperature, and negatively related to relative humidity and day of year.
Capture probability was lower if there was precipitation in the 24 h before trap-checking than if there was no precipitation (Figure 4).
Mean daily recapture probability was higher at Site I than at Site    Table 3). In contrast, the average probability of being available for capture was equal between the sexes at Site I (Pr[γ m,C > γ f,C ] = 0.45; Table 3).

| DISCUSS ION
Integrating models of distinct but related demographic processes can help extract more information from CMR data that are laborintensive to collect. Growth and survival are fundamental processes that interact to determine the age and size structure of a population, and therefore are natural choices to integrate into one model.
For our focal species, jointly analyzing growth and survival provided clear benefits. Given the low recapture rates of San Francisco gartersnakes, we could not measure the SVL of most individuals in the majority of years following their initial capture. Still, based on the importance of size to the demography of vertebrates (Sauer & Slade, 1987a) and for survival in snakes in particular (Brown & Weatherhead, 1999;Hyslop et al., 2012;, it was desirable to account for snake size when modeling survival.
Embedding a growth model into our CMR model allowed us to impute snake SVL when it was not observed, while ensuring that uncertainty in this state variable was propagated when estimating relationships between size and survival.
Although there was substantial uncertainty in annual survival for the smallest and largest snakes (due to small sample sizes and less precise estimates of availability for capture), there was a pattern of nearly equal survival for snakes between 300 and 550 mm SVL, and a decrease in survival for female snakes as SVL increased from 550 mm to 700 mm. Whether survival further decreased or plateaued for the largest San Francisco gartersnakes was unclear due to our small sample of snakes >700 mm SVL. Reduced survival for larger female San Francisco gartersnakes could be a sign of senescence (Reinke et al., 2020) or the survival costs of reproduction for large females (Luiselli et al., 1996;Madsen & Shine, 1993). Previous studies have documented costs of reproduction in female natricine snakes such as reduced locomotor ability (Seigel et al., 1987) and TA B L E 3 Parameters, prior distributions, and summary of posterior distributions from the multistate robust-design model for San Francisco gartersnakes (Thamnophis sirtalis tetrataenia) captured at two sites from 2007 to 2020. The prior "N(mean, SD)" represents a normal distribution with mean and standard deviation (SD), "Exp(scale)" is an exponential distribution with a scale parameter, a "Beta(shape1, shape2)" is a beta distribution with two shape parameters, and Uniform(lower, upper) is a uniform distribution between a lower and upper limit. Mean and SD are summary statistics from the posterior distribution for each parameter, and 2.50% and 97.50% indicate percentiles of the posterior distribution lower post-partum body condition (Brown & Weatherhead, 1997).
Actuarial senescence (a decrease in survival with age) has been documented for colubrid snakes  including some populations of the western terrestrial gartersnake (T. elegans; Colchero et al., 2019;Miller et al., 2011). The strength of actuarial senescence can depend on where a population falls on the fast-slow continuum of life history, with greater senescence expected in populations with fast life histories .
Although the San Francisco gartersnake has a restricted range in coastal northern California, our results hint that populations might differ in how each sex balances the needs for growth and survival to sexual maturity. Within our study populations, there were a few differences in patterns of growth and survival of female and male snakes. The finding that females reached greater lengths than males in each population is in keeping with other studies of San Francisco gartersnakes (Barry, 1994;Reeder et al., 2015). Males reached larger asymptotic length at Site I than at Site C, but females had similar asymptotic lengths at all five sites. Females exhibited higher annual survival than males at Site C, and males at Site C approached their asymptotic length at a faster rate than females. At Site I, the two sexes had nearly equal rates of growth, with some support for males having higher survival than females. We observed more spatial than temporal variation in snake growth rates, indicating that sites explain more of the variation in growth than annual fluctuations in the environment. In contrast, for survival, there were greater differences between the two sexes within a site than between mean estimates at each site. Data on reproductive traits from multiple populations could help elucidate the degree of life-history variation in San Francisco gartersnakes. It is possible that differences in climate between coastal and inland populations could be one factor in determining life-history strategies of San Francisco gartersnake populations. Daily temperature maxima and the influence of marine layer-induced fog can differ between inland and coastal sites on the San Francisco Peninsula (Torregrosa et al., 2016). Prey availability has also been identified as a cause of life-history variation among populations of the western terrestrial gartersnake (Bronikowski & Arnold, 1999), but we lacked quantitative data on the abundance and composition of prey at our sites.
Gartersnakes have been the subject of much research into lifehistory strategies and trade-offs between reproduction, growth, and survival. Substantial variation in life-history traits existed among populations of the common gartersnake (T. sirtalis) in Canada, with females from eastern populations reaching sexual maturity at a smaller size and producing larger litters of smaller neonates compared with western populations (Gregory & Larsen, 1993 The flexibility of using splines in a CMR model comes at a cost of needing to estimate additional parameters (e.g., compared with an intercept and slope in a linear function), and potentially greater uncertainty in the response variable where data are sparse. Overall, our recapture rates of San Francisco gartersnakes were low, with only 18% of snakes recaptured in at least one year following their initial capture and marking. This led to high uncertainty in our survival estimates, particularly for the youngest, smallest snakes and the largest snakes, which were rarely recaptured. For elusive snake species, complex statistical models alone cannot resolve the issue of low capture rates (Steen, 2010). Increasing recapture probabilities compared with the current study could require increased spatial coverage of trap arrays within a site (to reduce the effect of temporary emigration) and a longer duration of sampling within each year. For example, Reeder et al.
(2015) sampled a population of San Francisco gartersnakes with a much greater density of trap arrays and for a longer duration each year (>70 days) than in the current study, and were able to estimate abundance of snakes >300 mm SVL with good precision. Still, increased sampling effort (i.e., increased trap density) alone cannot overcome low recapture rates for individuals of certain life stages if the sampling method exhibits bias in the sizes captured Willson et al., 2008), or if the sampling period is not aligned seasonally F I G U R E 7 Estimated relationship between snake snout-vent length (SVL; mm) and availability for capture (γ) of San Francisco gartersnakes (Thamnophis sirtalis tetrataenia) from the robust-design model. Panel a presents estimated relationships from Site C, and panel b presents estimated relationships from Site I. Solid lines are posterior means, and shaded regions are 95% credible intervals. Black lines and gray shading represent female snakes, and red lines and shading represent male snakes. Vertical tick marks on the x-axis represent the SVL of captured snakes of each sex, shifted slightly to improve visibility with the availability of certain life stages (e.g., if sampling takes place before parturition and therefore neonates are unavailable).
Reducing uncertainty in our estimates of survival for the smallest and largest San Francisco gartersnakes could require a combination of increased sampling effort, greater duration of sampling within a year, and complementary approaches to studying the demography of these snakes. For reptiles, early life stages are notoriously difficult to capture and recapture (Pike et al., 2008). Estimating survival rates for small, young snakes could require integrating fecundity data along with CMR into an integrated population model (Schaub & Abadi, 2011). Likewise, estimates of availability of San Francisco gartersnakes for capture exhibited large uncertainty. The use of spatial capture-recapture models that account for snake movement (Ergon & Gardner, 2014) could improve estimates of availability for capture and lead to more precise estimates of survival. Unfortunately, within-year captures of San Francisco gartersnakes in our data were too sparse to fit spatial capture-recapture models and estimate the locations of individuals' activity centers. With increased sampling effort at a site, it might be possible to obtain enough within-year recaptures to fit spatial capture-recapture models, at least for snakes in the size range best sampled with funnel traps. Another potentially fruitful approach might be to use radiotelemetry of large adult females, in conjunction with survival analysis (Williams et al., 2002), to precisely estimate annual survival rates of this life stage.

| Summary and conclusions
Our study is an important step for developing a robust predictive framework for San Francisco gartersnake population dynamics.
Our results have shown that individual size influences growth and survival rates in this endangered subspecies and that growth and survival vary among populations and over time. The growth and survival parameters estimated here, in conjunction with fecundity data, could inform a size-structured demographic model that would identify size ranges of San Francisco gartersnakes to focus conservation efforts on to increase population viability. An important next step would be to link observed individual, interpopulation, and temporal variation in growth and survival to environmental variables and land management practices. Researchers with long-term CMR data could benefit from integrating models of growth, survival, and other demographic processes to reveal insights into how individual covariates affect vital rates and ultimately population growth.

ACK N OWLED G M ENTS
This work was funded by the U.S. Fish and Wildlife Service, the U.S.
Geological Survey Ecosystems Mission Area, and the Midpeninsula